library(ggplot2)
require(data.table)
library(ncdf4)

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"

crops = c("Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean")
case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")



#===============Regression analysis of CO2 effect==================
yield6crop.all <- readRDS(paste0(wd1,"/data/RData/yield6crop.all.Rds")) #from "Fig.1e.Global_yield_trend-final.R"

df.yield <- copy(yield6crop.all)
df.yield <- df.yield[, .(Prod=sum(Yield..MtC.year.), Area=sum(Area..Mha.)), by =.(Case,Crop,Year)]
df.yield <- df.yield[, Yield:=Prod/Area/0.45*0.85] #convert yield to tonnes/ha 
unique(df.yield$Crop)
unique(df.yield$Year)
cropname = unique(df.yield$Crop)


#CO2 difference between rcp45 and rcp 85
fco2_rcp45 <- nc_open(paste0(wd1,"/data/NCL/fco2_datm_rcp4.5_1765-2500_c130312.nc"))
fco2_rcp85 <- nc_open(paste0(wd1,"/data/NCL/fco2_datm_rcp8.5_1765-2500_c110919.nc"))
CO2_45=ncvar_get(fco2_rcp45, "CO2")[241:334] #2006-2099
CO2_85=ncvar_get(fco2_rcp85, "CO2")[241:334] #ppmv

co2.ts = data.table(co2= c(CO2_45, rep(CO2_85,times=4)), Year=rep(2006:2099,times=5), 
                    Case=rep(c("RCP45_85LU","RCP85","RCP85_SAI","RCP85_MSB","RCP85_CCT"), each=94)) 


#==============Step1: isolate CO2 fertilization effect

#use absolute CO2 change, easier to remove CO2 effect later no matter which case is the reference
#================co2 fertilization effect from low to high
#(base case RCP45, from RCP45 to RCP85)
co2.dif1 = (CO2_85 - CO2_45) #relative to dynamic CO2 of RCP45

co2.tr45 = (CO2_45 - CO2_45[1]) #CO2 trend of RCP45s relative to 2006
co2.tr85 = (CO2_85 - CO2_85[1]) #CO2 trend of RCP85s relative to 2006
co2.dif2= co2.tr85-co2.tr45 #use difference in co2 trends relative to 2006 level when detrending co2 effect from all cases
co2.dif2= co2.dif1 #identical as above 

co2.R2=co2.slope=co2.sd=co2.slope1=co2.slope2=co2.sd1=co2.sd2=co2.R22=dif.all=c()
for (ncft in 1:6) {
  yield.dif = df.yield[Crop==cropname[ncft] & Case=="RCP85",]$Yield - df.yield[Crop==cropname[ncft] & Case=="RCP85_45CO2",]$Yield
  #because land use changes every year in both RCP45 and RCP85, has to use % change relative to timeseries of RCP85_45CO2 (base)
  pct.dif = yield.dif/df.yield[Crop==cropname[ncft] & Case=="RCP85_45CO2",]$Yield*100 
  #the CO2 fertilization effect (from low CO2 to high CO2) will be halved from all cases for soy
  dif.all <- c(dif.all, pct.dif) #better use AvgYield to predict effect so that coefficients do not depend on the ref case
  #df = data.frame(yield.dif, co2.dif=CO2_45 - CO2_85)
  #mod.co2 <- lm(yield.dif ~ co2.dif, data=df) #absolute effect on yield
  df = data.frame(pct.dif, co2.dif = co2.dif1) #relative to dynamic CO2 of RCP45
  lm.co2 <- lm(pct.dif ~ co2.dif + 0, data=df) #no intercept (co2 is the only change factor between RCP85_45CO2 and RCP85)
  co2.slope <- c(co2.slope, lm.co2$coefficients[1]) #mean effect
  co2.sd <- c(co2.sd, summary(lm.co2)$coefficients[1,2]) #Std.Error
  co2.R2 <- c(co2.R2, summary(lm.co2)$r.squared)
  
  df2 = data.frame(pct.dif, co2.dif=co2.dif2, co2.tr45, co2.tr85) #relative to begining level of RCP45/85, same as co2.dif1
  #lm2.co2 <- lm(pct.dif ~ co2.dif +I(co2.dif^2) +0 , data=df2) #pct effect: use quadratic term as all crops show acclimation to CO2 fertilization
  lm2.co2 <- lm(pct.dif ~ co2.dif+I(co2.tr85^2-co2.tr45^2) +0 , data=df2)
  co2.slope1 <- c(co2.slope1, lm2.co2$coefficients[1])
  co2.slope2 <- c(co2.slope2, lm2.co2$coefficients[2])
  co2.sd1 <- c(co2.sd1, summary(lm2.co2)$coefficients[1,2])
  co2.sd2 <- c(co2.sd2, summary(lm2.co2)$coefficients[2,2])
  co2.R22 <- c(co2.R22, summary(lm2.co2)$r.squared)
}
#fertilization effect of co2 for each crop
cropname[1:6];co2.slope;co2.sd;co2.R2  #Soybean show saturation of CO2 fertilization, use quadratic term 
cropname[1:6];co2.slope1;co2.slope2;co2.R22 

summary(co2.dif1) #absolute max 380ppm, mean 108ppm; 
#co2 fertilization effect per 200 ppm (same as FACE experiment) increase of CO2
cropname[1:6];co2.slope*200 #(similar to FACE)
cropname[1:6];co2.tr85[73]-co2.tr45[73];co2.slope1*200+co2.slope2*(co2.tr85[73]^2-co2.tr45[73]^2)
#6.258913 9.893404 13.669582 15.544782 19.399228 44.595403% for maize, sugarcane, rice, cotton, wheat, and soy    

#co2 fertilization effect for doubled CO2 (380ppm, 70%)
cropname[1:6];co2.slope*380
#co2 fertilization effect if considering change relative to 2006 level
co2.slope*co2.tr85[94] #soybean > 100%
cropname[1:6];co2.tr85[94]-co2.tr45[94];co2.slope1*380+co2.slope2*(co2.tr85[94]^2-co2.tr45[94]^2)
#9.16835 16.45054% for maize and sugarcane, 30.67034 25.43299 20.15886% for wheat, cotton, rice, 74.79507% for soy

summary(co2.tr85);summary(co2.tr45) #pct change relative to 2006 level
summary(co2.dif2) #pct difference relative to 2006 level: mean 28%, max 100%
#soybean has too high CO2 fertilization effect
co2.dif2*co2.slope1[3] + (co2.tr85^2-co2.tr45^2)*co2.slope2[3]
#co2.slope1[3]*380+co2.slope2[3]*380^2 #use co2.dif^2 
crops[3];co2.slope1[3]*380+co2.slope2[3]*(co2.tr85[94]^2-co2.tr45[94]^2) #use co2.tr85^2-co2.tr45^2; very similar as co2.dif^2



#==========bias correction for soybean; (better to use AvgYield to predict and remove CO2 effect from all different cases)
#halving the overly high CO2 fertilization effect (e.g. Method 1: rcp45.yield/(1+co2.tr45.corr/100*0.5))
co2.dif.corr <- co2.slope1[3]*co2.dif1 + co2.slope2[3]*((co2.tr85)^2 - (co2.tr45)^2)
co2.tr45.corr <- co2.slope1[3]*co2.tr45 + co2.slope2[3]*(co2.tr45)^2
co2.tr85.corr <- co2.slope1[3]*co2.tr85 + co2.slope2[3]*(co2.tr85)^2
mean(co2.dif.corr[85:94]);mean(co2.tr45.corr[85:94]);mean(co2.tr85.corr[85:94])

#============End co2 correction for soy




#==========Plot original and corrected CO2 effect on average yield (quadratic function)
df.co2 = data.table(co2.effect.m1=rep(co2.slope1, each=94), co2.effect.sd1=rep(co2.sd1, each=94), #effect per 1ppm CO2 increase
                    co2.effect.m2=rep(co2.slope2, each=94), co2.effect.sd2=rep(co2.sd2, each=94),
                    R2 = rep(co2.R2, each=94),
                    Crop=rep(cropname[1:6], each=94), Year=rep(2006:2099, 6)) #6 crops x 94 years; 
df.co2$co2.dif=rep(co2.dif2,6) #modeled dif
df.co2$co2.tr=rep(seq(0, 400,length=94),6)  #predict 
df.co2$co2.tr45=rep(co2.tr45,6) 
df.co2$co2.tr85=rep(co2.tr85,6) 
df.co2$yield.dif=dif.all
df.co2 <- df.co2[,effect:=co2.dif*co2.effect.m1 + (co2.tr85^2-co2.tr45^2)*co2.effect.m2]

soy_cor <- copy(df.co2[Crop=="Soy"])
soy_cor[Crop=="Soy",yield.dif:=yield.dif-(co2.dif*co2.effect.m1 + (co2.tr85^2-co2.tr45^2)*co2.effect.m2)*0.5] #bias correction2
soy_cor[Crop=="Soy",effect:=effect/2]
soy_cor[Crop=="Soy", Crop:="Soy_halved"]
df.co2_cor <- rbind(df.co2, soy_cor)

save(list = c("df.co2","df.co2_cor"), file = paste0(wd1,"/data/RData/Fig.ED3.CO2_effect.RData")) #statistical prediction matches well with modeled difference













